#ifndef ESDIRK_HPP
#define ESDIRK_HPP

#include "IRK.hpp"

class ESDIRK : public IRK {
public:
  ESDIRK(int dim, double dt, double t0, double t1) : IRK(dim, dt, t0, t1) {
    m_s = 6;
    initializeSolver();
  }

  bool initializeSolver() {
    m_A.resize(6, 6);
    m_b.resize(6);
    m_b << 82889.0 / 524892.0, 0, 15625.0 / 83664.0, 69875.0 / 102672.0,
        -2260 / 8211.0, 0.25;
    m_c.resize(6);
    m_c << 0, 0.5, 0.5, 1, 1, 1;

    // Initialize all elements to 0 first
    m_A.setZero();

    // Fill in the non-zero elements
    m_A(1, 0) = 1.0 / 4.0;
    m_A(1, 1) = 1.0 / 4.0;

    m_A(2, 0) = 8611.0 / 62500.0;
    m_A(2, 1) = -1743.0 / 31250.0;
    m_A(2, 2) = 1.0 / 4.0;

    m_A(3, 0) = 5012029.0 / 34652500.0;
    m_A(3, 1) = -654441.0 / 2922500.0;
    m_A(3, 2) = 174375.0 / 388108.0;
    m_A(3, 3) = 1.0 / 4.0;

    m_A(4, 0) = 15267082809.0 / 155376265600.0;
    m_A(4, 1) = -71443401.0 / 120774400.0;
    m_A(4, 2) = 730878875.0 / 902184768.0;
    m_A(4, 3) = 2285395.0 / 8070912.0;
    m_A(4, 4) = 1.0 / 4.0;

    m_A(5, 0) = 82889.0 / 524892.0;
    m_A(5, 2) = 15625.0 / 83664.0;
    m_A(5, 3) = 69875.0 / 102672.0;
    m_A(5, 4) = -2260.0 / 8211.0;
    m_A(5, 5) = 1.0 / 4.0;

    return true;
  }
};

#endif